EXACT AND NON-STIFF SAMPLING OF HIGHLY OSCILLATORY SYSTEMS: 
AN IMPLICIT MASS-MATRIX PENALIZATION APPROACH 
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Abstract. We propose and analyze an implicit mass-matrix penalization (IMMP) technique which enables efficient 
and exact sampling of the (Boltzmann/Gibbs) canonical distribution associated to Hamiltonian systems with fast degrees of 
freedom (fDOFs). The penalty parameters enable arbitrary tuning of the timescale for the selected fDOFs, and the method 
is interpreted as an interpolation between the exact Hamiltonian dynamics and the dynamics with infinitely slow fDOFs 
(equivalent to geometrically corrected rigid constraints). This property translates in the associated numerical methods into a 
tunable trade-off between stability and dynamical modification. The penalization is based on an extended Hamiltonian with 
artificial constraints associated with each fDOF. By construction, the resulting dynamics is statistically exact with respect to 
the canonical distribution in position variables. 

The algorithms can be easily implemented with standard geometric integrators with algebraic constraints given by the 
expected fDOFs, and has no additional complexity in terms of enforcing the constraint and force evaluations. The method is 
demonstrated on a high dimensional system with non-convex interactions. Prescribing the macroscopic dynamical timescale, 
it is shown that the IMMP method increases the time-step stability region with a gain that grows linearly with the size 
of the system. The latter property, as well as consistency of the macroscopic dynamics of the IMMP method is proved 
rigorously for linear interactions. Finally, when a large stiffness parameter is introduced, the IMMP method can be tuned to 
be asymptotically stable, converging towards the heuristically expected Markovian effective dynamics on the slow manifold. 

Key words. Hamiltonian systems, NVT ensemble, slow/fast systems, Langevin dynamics, constrained dynamics, hybrid 
Monte Carlo. 
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1. Introduction. This paper deals with stochastic numerical integration and sampling of Hamil- 
tonian systems with multiple timescales. The main motivation is to develop numerical methods which 
sample exactly canonical distributions, while integrating the fastest degrees of freedom only statistically. 
Furthermore, one also seeks good approximation of dynamical behavior, at least at large temporal scales. 

Hamiltonian systems with multiple timescales typically appear in molecular dynamics (MD) simula- 
tions, which have become, with the aid of increasing computational power, a standard tool in many fields 
of physics, chemistry and biology. However, extending the simulations to physically relevant time-scales 
remains a major challenge for various large molecular systems. Due to the complexity of implicit methods, 
the time scales reachable by standard numerical methods are usually limited by the rapid oscillations of 
some particular degrees of freedom. Since the sampling dynamics has to be integrated for long times, the 
time-step restriction associated with fast oscillations/short time scales in molecular systems contributes to 
the high computational cost of such methods. Yet the physical necessity of resolving the fast degrees of 
freedom in simulations is often ambiguous, and efficient treatment of the fast time scales has motivated 
new interest in developing numerical schemes for the integration of such stiff systems. 

The problem of integrating stiff forces is relevant both for the direct numerical simulation of the 
Hamiltonian dynamics, and for the less restrictive problem of designing a process that samples the canon- 
ical ensemble. Sampling from the canonical distribution can be achieved by Markov chain Monte Carlo 
(MCMC) algorithms based on an a priori knowledge of possible moves, with a Metropolis-Hastings ac- 
ceptance/rejection corrector (a historical reference is |30j). For complex molecular systems however, such 
global moves remain unknown in general, and sampling methods consists generically in using either an 
Hamiltonian dynamics with a thermostat (e.g., a Langevin process), or a drifted random walk (Brownian 
dynamics), corresponding to over-damped Hamiltonian dynamics (see for a review and references on 
classical sampling methods). Brownian dynamics of systems with multiple time scales suffers from similar 
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stability restrictions (see |201 139] for some practical issues on Brownian dynamics simulation for Molecular 
Dynamics) . 

Broadly speaking, one may start by recognizing two approaches to the numerical treatment of stiff 
systems: 

(i) Semi-implicit, multi-step integrators and their variants (e.g., the textbooks [311 E], or the review paper 

[lOj and references therein, [20j for Brownian dynamics), which attempt to resolve microscopic 
highly oscillatory dynamical behavior. 

(ii) Methods with direct constraints, where the highly oscillatory degrees of freedom are constrained to 

their equilibrium value (e.g., [Ml [251 ES] and references therein). 
In spite of their differences the common key feature of all these methods is to balance a trade-off between 
stability restrictions and implicit time-stepping form, or in other words, between the computational effort 
associated with small time steps, and the computational cost of solving implicit equations implied by the 
stiffness. 

The method developed in this paper aims at giving an appropriate interpolation between exact dy- 
namics and constrained dynamics considered in the second family of methods mentioned above. Although 
constrained dynamics remove, in principle, the stiffness of the associated numerical scheme, it introduces 
new difficulties and numerical problems. As an approximation to the original dynamics it modifies im- 
portant features of the system; most importantly, the original statistical distribution. The principal goal 
of the proposed method is to replace direct constraints by implicit mass-matrix penalization (IMMP), de- 
tailed in Section [3l which integrates fDOFs, but with a tunable mass penalty. This approach guarantees 
that the canonical distribution is computed exactly, and that a freely tunable trade-off between dynamical 
modification and stability can be obtained. 

The idea of adjusting mass tensors in order to slow down fast degrees of freedom goes back to 3 . In this 
paper, the author proposes to modify the mass tensor with respect to the Hessian of the potential energy 
function in order to confine the frequency spectrum to low frequencies only. Two natural drawbacks of this 
procedure arise from the costly computation of the second-order derivatives of the potential, and from the 
bias introduced when the adjusted mass-tensor is adapted during the dynamics. Such an approach seems 
inevitable when the fDOFS are unknown, but in many cases, the fast degrees of freedom are explicitly given 
by the structure of the system (e.g., co-valent and angle bonds in molecular chains). To our knowledge, 
mass tensor modification have been used in practical MD simulations by increasing the mass of some well- 
chosen (e.g., light) atoms |261 129j . The aim of this paper is to propose a more systematic mass-tensor 
modification strategy. 

The proposed method relies on the assumption that the system Hamiltonian is separable with quadratic 
kinetic energy 



and that the "fast" degrees of freedom (fi,..,Cn) are explicitly defined smooth functions of the system 
position 



We emphasize that the knowledge of "fast forces" is not required, and the variables ^ can be chosen 
arbitrarily. If the fDOFs are not identified the method retains its approximation properties while not 
performing efficiently. The fDOFs are penalized with a mass-tensor modification given by 



Hip,q)^-p^M-'p + V{q) 



(1.1) 



q = (qi, ...,qd)^-^ (fi('?), ■■■,S.n{q)) ■ 



(1.2) 



(1.3) 



where v denotes the penalty intensity, and a "virtual" mass matrix associated with the fDOFs. The 
modification does not impact motions orthogonal to the fDOFs. The position dependence of the mass- 
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penalization introduces a geometric bias. This bias is corrected by introducing an effective potential 



which will turn out to be a i-perturbation of the usual Fixman corrector (see |19j ) associated with the 
sub-manifolds defined by constraining the fDOFs f . The key point is then to use an implicit representation 
of the mass penalty with the aid of the extended Hamiltonian 



The auxiliary degrees of freedom z are endowed with the "virtual" mass- matrix M^. The constraints (Ci/) 
are applied in order to identify the auxiliary variables and the fDOFs ^ with a coupling intensity tuned 
by V. The typical time scale of the fDOFs is thus enforced by the penalty v. The system is coupled to 
a thermostat through a Langevin equation (|2.3p . which yields a stochastically perturbed dynamics that 
samples the equilibrium canonical distribution. We then obtain the following desirable properties: 

1. The associated canonical equilibrium distribution in position is independent of the penalty v. 

2. The limit of vanishing penalization {y = 0) is the original full dynamics, enabling the construction 
of dynamically consistent numerical schemes. 

3. The limit of infinite penalization is a standard effective constrained dynamics on the "slow" man- 
ifold associated with stiff constraints on ^. 

4. Numerical integrators can be obtained through a simple modification of standard integrators for 
effective dynamics with constraints, with equivalent computational complexity. 

The dynamics associated with the IMMP Hamiltonian (|1.5p is detailed in (|3.6p . see Section [3l Its 
numerical discretization by a leap frog/ Verlet scheme with constraints (usually called "RATTLE") is given 
by (|4.2p . The flow can then be corrected by a Metropolis step to obtain exact canonical sampling, which is 
reminiscent of so-called Hybrid Monte-Carlo methods (see references in Section U) . When this correction 
is introduced, the gradient of the Fixman potential (|1.4p need not be computed. The numerical aspects of 
the method are presented in Section [H 

By including a penalty, the proposed method modifies the original Hamiltonian. However, the mass 
penalty can be also thought of as depending on the time step v := v{8t") leading to consistent schemes. In 
Section [S] and in Section [71 the penalty (and the fastest timescales) will even grow to infinity v — > -foo as 
the dimension of the system — > cx) or a stiffness parameter e ^ 0. In both cases, we encounter generic 
situations where the IMMP method leads to dynamically consistent limits, a partial differential equation, 
and an effective dynamics on a "slow" manifold, respectively. 

High-dimensional systems usually contain a large variety of timescales, and are therefore challenging 
test cases. A linear chain of particles with repulsive non-convex interactions is numerically studied in 
Section [5l It is shown that a large mass-penalty of order v :— DN, where A^ is the system size, induces a 
gain in time stepping of order A^, while numerical evidence are given that the macroscopic dynamics, given 
by a formal stochastic partial differential equation, remains of order 1 (in particular, the convergence to 
equilibrium). Rigorous proofs with explicit scalings of these behaviors is provided for the linear case (i.e. 
harmonic interactions) in Section [6l and consistence of the IMMP macroscopic dynamics when the scaled 
penalty vanishes — s- is even demonstrated. 

In Section [7l we introduce the stiffness parameter e and show that the penalty intensity can be scaled 
as v :— V 1 1 in order to obtain asymptotically stable dynamics in the limit e ^ 0. We prove that the 
dynamics converge towards the expected Markovian effective dynamics on the slow manifold. We also 
present analysis of stability properties of the proposed scheme. 
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V^fix..(g) = — ln(det(Af,(g))) , 
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2. Formulation. We consider a Hamiltonian system in the phase-space M'' x R'^ with the Hamiltonian 
H in the form 

H{p,q)^^p^M-'p + V{q), (2.1) 

We use generic matrix notation, for instance the Euclidean scalar product of two vectors pi,P2 G de- 
noted by pjp2, and the gradients of mappings from to K" with respect to standard bases are represented 
by matrices 

i^qOij = i^gOji = 1^ , i = 1, . . j = 1, . . . 

When the system is thermostatted, i.e., kept at the constant temperature, the long time distribution 
of the system in the phase-space is given by the canonical equilibrium measure at the inverse temperature 
/3 (also called the NVT distribution) given by 

fj,{dpdq) = ^e~'^"^P^'^Updq, Z^f e-^^^P^^Updq, (2.2) 

with the normalization constant Z < oo. The standard dynamics used to model thermostatted systems 
are given by Langevin processes. 

Definition 2.1 (Langevin process). A Langevin process at the inverse temperature (3 with the Hamil- 
tonian H{q,p), {p,q) E M.'^ X R'', the d x d dissipation matrix 7, and the fluctuation matrix a is given by 
the stochastic differential equations 

( q ^ WpH 

{ ■ (2.3) 

[ p = -VgH --fq + aW , 

where is a standard white noise (Wiener process), and a G M'' x satisfies the fluctuation- dissipation 
identity 

For any 7, the process is reversible with respect to the stationary canonical distribution h2.}3(l . Furthermore, 
if 7 is strictly positive definite, the process is ergodic. Throughout this paper, the usual global Lipschitz 
conditions (see [31]) on H and ^ are assumed, ensuring well-posedness of the considered stochastic differen- 
tial equations. The analysis presented in the paper can be generalized to a position dependent dissipation 
matrix 7 = 7(9). 

The mapping ^ : M'' M", defines n < d degrees of freedom, given by smooth functions taking values 
in a neighborhood of 0. We assume that the mapping ^ is regular (i.e., with a non-degenerate Jacobian) 
in an open ^-neighborhood Og = {q \ 11^(9)11 < ^} of C~^(0)j hence defining a smooth sub-manifolds of K'^ 
denoted M.z = for z in a neighborhood of 0. The dependence of the potential with respect to the 

degrees of freedom ^ is expected to be "stiff" in the second variable. In Section [7] we will introduce the 



^Throughout the paper, stochastic integrands have finite variation, so that stochastic integration (e.g. Ito or Stratonovitch) 
need not be specified. 



Sampling highly oscillatory systems. 



5 



stiffness parameter e. We will assume in this section that such parameter dependence can be explicitly 
identified, and that the potential energy V can be written in the form 



where the function [/ : x M" M satisfies the coercivity condition lim^^+oc U{q,z) — +00. The fast 
degrees of freedom ^ of states at a given energy then remain in a closed neighborhood of the origin as the 
stiffness parameter e — > 0. In this limit the system is confined to the sub- manifold A4q, which is usually 
called the "slow manifold" . 

3. The implicit mass-matrix penalization method. In this section we focus on properties of 
the IMMP method. The multiscale structure of the potential V need not be known in order to apply 
the method. Thus, in this section, we consider the potential V in the form where we do not impose the 
structural assumption (|2.4p on the potential function y : K'' ^ R. 

The new, penalized mass- matrix of the system is the position dependent tensor defined in (jl.3p . The 
associated modified impulses are denoted 



When v becomes large, the velocities are bound to remain tangent to the manifolds {q\^{q) ~ z}, and 
orthogonal motions are arbitrarily slown down. Conversely, when v — 0, one recovers the original highly 
oscillatory system. Since the modification in AIi, depends on the position q, new geometry is introduced 
and an additional correction (|1.4p in the potential energy is required in order to preserve original statistics 
in the position variable. This correction is in fact close to the standard Fixman corrector for v large (see 
(|7.5p ). Defining G{q) as the n x n Gram matrix associated with the fast degrees of freedom 




(2.4) 



= M^{q)M-^p. 



(3.1) 



G(g)=V^^A/-i%e, 



(3.2) 



one has the following property of the correcting potential. 
Proposition 3.1. Up to an additive constant, we have 




(3.3) 



and thus (up to additive constants) 




lim Pfix,!/ — Vfix — Indet {G{q)) , and lim Va^.i, = . 



Proof. Using the identity for a non-diagonal matrix J of dimension ni x 712: 



det (Id„i + JJ^) = det (Id„2 4- J) , 



one observes 



det (AU) = det (M) det (i^^M^ det (G + ^AC^) 



from which the expression for the corrected Fixman potential follows. □ 
The associated modified Hamiltonian is then given by 



H^ip^, q) = -pIAU ^Pu + V{q) + Vfix,^(g) , 



(3.4) 
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and Hq — H is the original Hamiltonian p.ip . 

By construction, statistics of positions q of the mass penalized Hamiltonian are independent of the 
penalization, leading to the exact canonical statistics in position variables. 

Proposition 3.2 (Exact statistics). The canonical distribution associated with the mass-penalized 
Hamiltonian p.4p is given by 

H^idp, dq) = ^e-''^-'(P-9)rfp^ dq . (3.5) 
Its marginal probability distribution in q is 

1 f ^-BH.Jr,.,.^).^^ e-'^VWdg 



i_ J e-f^"^(P-'^Up, 



Je-PV(q) dq 

which is the original canonical distribution \2.2\ in the position variables, and is independent of the mass 
penalization parameter v. 

Proof. The normalization of Gaussian integrals in the py variables yields 

which is cancelled out by the Fixman corrector Vrx,!/ and the result follows. □ 

Sampling such a system can be done using the standard Langevin stochastic perturbation as detailed 
in Definition 12.11 However, the direct discretization of the equation of motion given by (e.g., by an 
explicit scheme) is bound to be unstable (from non-linear instabilities) when the fast degrees of freedom are 
not affine functions. In order to construct stable schemes one may rather use an implicit formulation of the 
Hamiltonian p.4p . in conjunction with a solver which enforces the constraints. To obtain such a formulation 
we extend the state space with n new variables (zi,..,z„), and associated moments (pz-^, ..,pz„). The 
auxiliary mass-matrix for the new degrees of freedom is then given by ■ The new extended Hamiltonian 
of the system i/iMMP, defined by (|1.5[) . is now defined in R''+" x R''+", where n position constraints denoted 
by (Cjy) are included. This construction implies n hidden constraints on momenta. The equivalence of the 
two Hamiltonians (j3.4p and (jl.Sp formulations is stated as a simple separate lemma. 

Lemma 3.3. The equations of motion associated with the penalized mass-matrix Hamiltonian p.4p or 
the extended Hamiltonian with constraints (jl.Sp are identical. 

Proof. The Lagrangian associated with Himmp is given by 

LiMMp{q, z, q, z) = ^q^Mq + ^z'^M^z - V{q) - Vfi^,t,{q) , 

and includes hidden constraints on velocities i = v'^qS.q implied by the constraints {C^) on position 
variables. Replacing i and z in Ximmp by their expressions as functions of q and g, one obtains the 
Lagrangian associated with H^. □ 

The stochastically perturbed equations of motion of the Langevin type associated with (|1.5p define the 
dynamics with implicit mass-matrix penalization. 

Definition 3.4 (IMMP). The implicit Langevin process associated with Hamiltonian Himmp and 
constraints {C^) is defined by the following equations of motion 

q = Af-ip 

z ^ M-^pz 

p = -VqV{q) - V,Ffi,,,((z) ^^q + aW- Vg^X 

^ (3.6) 

Pz = -IzZ -\- a-zWz -\- - 

V 

aq)--, (c.) 

V 
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The process W (resp. Wz ) is a standard multi- dimensional white noise, 7 (resp. j^) a dx d (resp. nx n) 
non-negative symmetric dissipation matrix, a (resp. a^) is the fluctuation matrix satisfying aa'^ — ^7 
(resp. azcr"^ — ^Iz)- The processes A G M" are Lagrange multipliers associated with the constraints (C^) 
and adapted with the white noise. 

This process is naturally equivalent to the explicit mass-penalized Langevin process in M."^ x M'' as- 
sociated with M,y. Moreover, when the penalization vanishes {ly 0), the evolution law of the process 
{Pt,(!t}t>o or {{pi^)t, <lt}t>o converges towards the original dynamics. 

Proposition 3.5. The stochastic process with constraints (j3.6p is well-posed and equivalent to the 
Langevin diffusion (see Definition \2.1\) . with the mass-penalized Hamiltonian H^, (j3.4p . and the 

dissipation matrix given by 

7.(9) =7 + '^'V,e7,V^C- 

Furthermore, the process is reversible and ergodic with respect to the canonical distribution (13. 5p (with 
marginal in position variables given by the original potential, i.e., up to the normalization, e~^^^'^^dq) . 
Proof. Imposing the constraints implies V^^il/~^p = ^M~^pz. Thus by the definition of p^ we have 

p„ ^p + vVqC pz . 

Since the position process {qt}t>o is of finite variation, a short computation shows that for each coordinate 
i = 1, ..,d 

.2 -T 



Furthermore, 



and thus 



Pl + :y{d,^y ^pz + v'q' y 18,^0 Pz ■ (3.7) 



Substituting the expressions for p and Pz from (j3.6p into p.7p we obtain 

Pv = -ip^V,M- V - VgF(g) - ^qV^-^A'i) - 79 - v^^qi IzZ + aW + vVqS, GzWz , (3.8) 
which yields the result. □ 

Proposition 3.6 (Small penalty). When v ^ the evolution law of the processes {pt,qt}t>o or 
{{Pi^)t,qt}t>a defined by the implicit equations (|3.6|) converges (in the sense of probability distributions on 
continuous paths endowed with the uniform convergence) towards the process solving the original Langevin 
dynamics (j2.3p . 

Proof. The stochastic differential equation defined hy q — M^^p^, and p.Sp has smooth coefhcients 
which depend on in a continuous fashion {v 1-^ M^, and v 1-^ Vflx,;^ are continuous). Standard results on 
weak convergence ([16]) of stochastic processes imply the result as stated. □ 

When the mass penalty tends to infinity, the IMMP process converges to a constrained process on the 
manifold Mzt=o — = zt=o} (surface measures are described in Section [X]). 

Proposition 3.7 (Large penalty). Consider a family of initial conditions indexed by v and satisfying 



sup \iy (C(9t=o) - zt=o)\ < +00 , 

V 
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and assume that the Gram matrix G is invertible in a neighborhood of A4zt=o- Then when v +00 the 
IMMP Langevin stochastic process p.6p converges weakly towards the decoupled limiting processes with 
constraints 



M- 



VgFfi^ -jq + crW- Vg^X , 



^{q) = zt=o , 
z = M-^p, , 
^ Pz = ~"fzz + azWz . 



(3.9) 



where {Xt}t>o o-f^ adapted stochastic processes defining the Lagrange multipliers associated with the con- 
straints (C). 

Furthermore, the process {qt,Pt\t>a defines an effective dynamics with constraints (see also Defini- 
tion^^ on the sub-manifold Aizt=o- reversible with respect to its stationary canonical distribution 
given, up to the normalization, by 



-P(H(p,q)+VnA<])) 



<^T-'M,^^g{dpdq) 



with the q-marginal e~^^^''^ S^{q)=o{dq). When 7 and 72 are strictly positive definite the process is ergodic. 

Proof. By a simple translation, it is sufficient to show the proposition for zt=o = 0. Satisfying 
the constraint (Ci/) in (|3.6p implies a hidden constraint in the momentum space, V^^Af^^p = ^M^^pz- 
Differentiating this expression with respect to time and replacing the result in p.6p yields an explicit 
formula for the Lagrange multipliers 



A 



with forces f^) 



(G+^M-^) 



1 



Hess (e) (M-ip, M~'p) + Vq^M-'f, - -M-'f, 



(3.10) 



fq = -V,V - %T4i,,, - jM-'p + aW , 
fz = -"fzM~'^pz + (JzWz, 

and the Hessian Hess {£,) of the mapping ^ acting on the velocities M~^p. This calculation shows that (|3.6p 
is in fact a standard stochastic differential equation with smooth coefficients, and thus has a unique strong 
solution. The coefficients of these stochastic differential equations are continuous in the limit ^ ^ 0, at 
least in a (S-neighborhood of Mq in which G is invertible. The formally computed limiting process is given 
by p.9p with the Lagrange multipliers solving 

A = (Hess (0 {M-^p, AT^p) + Vj^ M" V,) • 

By construction, this limiting process satisfies the constraint ^(q) = 0. Its coefficients are Lipschitz and 
the process is well-posed. As a result of those properties, the rigorous proof of weak convergence follows 
classical arguments, see |16j . that are divided into three steps 

(i) We truncate the process (|3.6p to a compact neighborhood oi Mo- 

(ii) The continuity of the Markov generator with respect to i implies tightness for the associated ^- 

scquences of truncated processes with the limit being uniquely defined by (|3.9p . 

(iii) The limiting process remains on A^o, which implies weak convergence of the sequence without trun- 

cation. 

The process (j3.9p is thus a Langevin process with constraints, exhibiting reversibility properties with respect 
to the associated Boltzmann canonical measure and is ergodic when 7 is strictly positive definite (see the 
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summary in Appendix [B| . Note that the g-marginal is geometrically corrected by the Fixman potential 
term. □ 

We conclude this section by discussing some consequences for numerical computations. 
Remark 3.8. 

(i) Proposition lS. 6\ and \3. 7| imply that IMMP scheme is a tunable interpolation between the exact dynamics 

()2.3p . and the dynamics with constraints (|3.9p . 

(ii) In Section^where a stiffness paramater, corresponding to a large penalty v, is introduced for slow/fast 

systems leads to a Markovian effective motion constrained on the slow manifold Mo, usually seen 
as a heuristically consistent Markovian approximation of the exact dynamics (see 133f ). 

4. Numerical integration. The key ingredient for achieving efhcient numerical simulation is to use 
an integrator that enforces the constraints associated with the implicit formulation of the mass penalized 
dynamics p.6p . The implicit structure of p.6p leads to numerical schemes that are potentially asymp- 
totically stable in stiff cases (Section [7]). On the other hand, when the penalization v vanishes with the 
time-step the scheme becomes consistent with respect to the original exact dynamics ()2.3p . One may then 
consider the mass-penalization introduced here as a special method of pre-conditioning for a stiff ODE 
system with an "implicit" , in the time evolution sense, structure. Here, the "implicit" structure amounts 
to solving the imposed constraints £,{q) = j; in p.6p . 

It lies outside the scope of this paper to make a review of standard numerical methods for constrained 
mechanical systems, we refer to 14 as a classical textbook, and to the series ( [38l [25l [HI ETJ [8] ) as a sample 
of works on practical developments of numerical methods. The IMMP method is presented with the classical 
leapfrog/ Verlet scheme that enforces constraints, usually called RATTLE in (|4.2p . It can be implemented 
by a simple modification of standard schemes constraining fDOFs. The scheme is second order, reversible 
and symplcctic. This choice is largely a presentation matter, for practical purposes one can refer to one's 
favorite numerical integrator for Hamiltonian systems with or without stochastic perturbations. 

As an option, one can also add a Metropolis acceptance/rejection corrector at each time step. If the 
underlying integrator is reversible and preserves the phase-space measure, this extension leads to a scheme 
which exactly preserves canonical distributions. Such algorithms are usually referred to as Hybrid Monte 
Carlo (HMC). Generically HMC are sampling algorithms which use the underlying dynamics of the system 
to generate moves in the configuration space, which are accepted or rejected according to a Metropolis rule. 
While this method may allow for larger time steps in the integrator it may also lead to higher computational 
costs when the acceptance ratio becomes low. The acceptance ratio often decreases when the dimension of 
the system increases However, many improvements and modifications ( [TTl [551 [51 [53] ) of the HMC 

algorithm have been developed since its introduction in [T51 [T3] • While the HMC methods were initially 
used in simulations applied to quantum statistical field theories, they have been also employed to wide 
range of simulations in macromolecular systems, e.g., 136] . 

We implement numerical discretization of the Langevin process with constraints (j3.6p obtained by 
splitting the Hamiltonian part and the Gaussian fluctuation/dissipation perturbation. The idea of using 
a Hybrid Monte Carlo step in discretizing a Langevin process goes back to |55]. Note that by using 
a Metropolis acceptance/rejection rule (HMC), or simply by weighting statistical averages, the forces 
associated with the Fixman corrector need not be computed when one is only interested in sampling. 

We recall that we consider the IMMP dynamics (|3.6p . which consists of a Langevin process with 
constraints defined by the following elements 

1. the Hamiltonian i?iMMP defined in (jl.Sp . 

■ fl o\ 

2. the dissipation matrix I ^ 'y ) ' 

3. the inverse temperature /3. 

Remark 4.1. A useful variant, when using HMC strategies, consists in modifying the Hamiltonian 
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^IMMP (|l-5p in the integrator (importance sampling) by neglecting the Fixman corrector 
HiMMp{p,Pz,q, z) = ip^M^V + ^pT^C + V{q) 

i{q) = - (a), 

where V is a potential that may he chosen arbitrarily (typically V = V). Indeed, only the underlying phase 
space structure, which does not depend on V , is necessary in HMC methods. The correct potential V^ + Vfix.jy 
has only to be used in the Metropolis step, or by weighting ensemble averages, to ensure exact canonical 
sampling. Thus potentially costly evaluations of the gradient of the Fixman corrector Vfix.i^ can be avoided. 

Scheme 4.1 (Leapfrog/Verlet algorithm for Langevin IMMP (|3.6p ). 
Step 1: Integrate the Hamiltonian part with: 

St 

Pn+l/2 = Pn ~ y + VgVfix,z.)(<J«) " ^q(,{q7i)K+l/2 

Pn+l/2 - Pn + -K+1/2 
qn+1 = qn + StM Pn+l/2 

Zn+1 = Zn + StM^ Pn+l/2 (4 2) 

' ^{qn+l) = ^ (^1/2) 

Pri+l = Pn+l/2 - y + VgVfix,,.) (<?„+! ) - Vg^(g„+i ) A„+i 

1 ^ 

Pn+1 = Pn+1 + - Wl 

V^^qn+i) A/- V+1 = ^M-i^+i (Ci) . 

Step 2: Integrate the Gaussian fluctuation/dissipation part with a mid-point Euler scheme with constraints 
(see Appendix\D^. 

To obtain exact sampling by correcting the time step errors, one needs to introduce a Metropolis correction. 
Scheme 4.2 (HMC). 

Step 1: Compute {q„+i, Zn+i,Pn+i,Pn+i) ™th and set AH^+i = HiuMpiqn+i, Zn+i,Pn+i,Pn+i) - 

^^IMMP {qn, ZmPmPn)- 

Step 2: Accept the step with the probability min(l, e~'^'^^"+^, otherwise reverse impulses and set 

{qn+l,Zn+l,Pn+l,Pn+l) = {qn, Zn, -Pn, -Pn) ■ 

Step 3: Integrate the Gaussian fluctuation /dissipation part with a mid-point Euler scheme (for details see 
Appendix\D\) . 

This numerically constructed Markov chain preserves the canonical distribution. 

Proposition 4.2 (Exact sampling). Assume that (14. 2p is defined globally. The numerical discretiza- 
tion of p.6p described in Scheme \4-.2\ generates a Markov chain that leaves the canonical distribution (|3.5p 
invariant, with marginal distribution in position variables is the original distribution e~^^^'^^dq, which is 
independent of the mass-penalization v . 

Proof. The statement follows from reversibility and measure preserving properties of Verlet schemes 
(see |14j)j from the Hybrid Monte Carlo rule (llSj), and from the construction of the mass-penalized 
Hamiltonian (Proposition [33]). □ 

One can now construct consistent schemes by letting the penalty v — vbt^ go to zero with the time-step, 
for some fc > 0. Indeed, Proposition 13.61 shows that the mass-penalized dynamics p.6p converges towards 
the exact original dynamics for v — ^ order . Consequently most of the usual numerical schemes will 
be consistent at their own approximation order, but bounded above by 2fc. Neglecting the order of the 
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fluctuation/dissipation part (Step 2 of Scheme [iTTj) . we deduce the following convergence property. The 
order of convergence refers to the maximal integer k such that the convergence of trajectories with respect 
to the uniform norm occurs at the rate 0{St^). 

Proposition 4.3 (Time-step consistency). Assume that the numerical flow (|4.2p is defined globally, 
and that v = vit^ . Then the IMMP numerical scheme (14. 2p is consistent of the order min(2, 2fc) with 
respect to the original exact limiting process (j2.3p without thermostat (7 = 0, a — 0). 

Proof. We consider the scheme (|4.2p in the variables {p, q, lypz, vz). Then the force field only depends 
on q and the global mass- matrix is smooth with respect to v^. The sub-manifold defining the constraints 
{i'^S,{q) = vz) is also smooth with respect to iP' . By the implicit function theorem we have that locally the 
RATTLE scheme is the standard leapfrog scheme (see and the local mapping depends smoothly 

on . Therefore the standard calculation of the order of the leapfrog scheme holds uniformly with respect 
to V, [Tlj. Now, the mass-penalized dynamics (|3.6p is a differential equation (deterministic here) whose 
coefRcients differ from the original process by a smooth perturbation of order . The result thus follows 
from applying a simple Gronwall argument. □ 

Note also that taking v +00 in (|4.2p gives a numerical scheme with constraints consistent with the 
limiting constrained dynamics of Proposition 13.71 For slow/fast systems (see Section [7^ . large penalty 
will lead to asymptotic stability in the stiffness limit. 

We conclude this section with a practical recipe for tuning the mass matrix penalty. This can be done 
for instance by computing the time averaged Metropolis acceptance ratio r = E [e"'^'-'^^"^] in Scheme l4?2l 
which gives a precise quantification of the time-step error. Then increasing the penalty v can save compu- 
tational time as long as it leads to an increase of the average ratio r. Indeed, this means that the selected 
fDOFS arc limiting the time-step stability region. Prescribing a time average ratio r (for instance 0.90), 
a maximal time-step dtmax associated with the largest penalty Vmax that is able to improve stability can 
be obtained in this way. Finally, one can set, for example, v = dt in Scheme 14.21 to obtain an order 2 
convergent scheme with increased stability region. 

5. A high dimensional numerical test case. The behavior of numerical methods for Hamiltonian 
systems becomes of particular interest when simulating high-dimensional systems. Such systems in general 
exhibit many different time-scales and one seeks to go beyond stability constraints implied by the fastest 
microscopic degrees of freedom. As mentioned in the introduction, globally implicit methods are usually 
too costly, while efficient splitting methods using minimal implicitness represent an active research area. 
In order to demonstrate applicability of the IMMP method we consider non-convex repulsive interactions 
between particles in the gas phase which makes implicit methods intractable. 

On the other hand, direct constraints on microscopic fast degrees of freedom are commonly used in 
practice and can introduce a potentially strong bias in the macroscopic behavior of the whole system. In the 
case presented here, a particle chain with pairwise interactions, one may be interested in some macroscopic 
observable, like the total length of the chain. Direct constraints would then lead to a totally rigid chain, 
losing completely the evolution of the latter macroscopic variable. 

We numerically study the IMMP strategy described in Section [3] and Section 31 and the latter is 
systematically compared with the exact dynamics numerically integrated by a simple Leapfrog/ Ver let 
scheme with a thermostat. Langevin thermostats are used, with a Metropolis/HMC step for the dynamics 
integrator as described by Scheme [4?2l An appropriate scaling of the chain of size N at the mass transport 
level with a large mass-penalty of order v :— vN are chosen. Numerical results show that IMMP induces a 
gain in time stepping of order N , while the macroscopic dynamics (in particular convergence to equilibrium), 
given by a formal non- linear stochastic partial differential equation, remains of order 1. Macroscopic 
dynamics of the IMMP method also stays close to the exact Verlet integration when re-scaled penalty v is 
small. 
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5.1. The particle chain model. The model we consider consists of a chain of particles which interact 
through a non-convex repulsive pairwise potential of the form 

1 Wi„t(r) = /(r) , r<ri„t, 

[ Wint(r) = , r > Tint , 

where f{rint) = /'(nnt) — and lim /(r) = +oo. Each particle is also individually submitted to 

r — ' — oo 

a macroscopic confining exterior potential. After converting to the non-dimensional form the typical 
quantities involved in the model enable us to write a scaling at the mass-transport level where the dynamics 
of the chain is described by the Hamiltonian 

^ N~l N 

HN{q.p) = ^p'^P + «mt(V- g) + woxt(ft) , (5.2) 

i=l i=l 

and by a coupling with an exterior thermal bath at the re-scaled inverse temperature 

(3n ■■= PN^^ . (5.3) 

In the expression (|5.2p the functions r G R i-^ Vint{r) E M and q h-> Vcxt{q) G R are the smooth interac- 
tion potential and the exterior potential rrespectively The linear operator : K^^^ having the 
components 

^ <l^+l - 1i -^^ 

^iH l/N ' ' -'-J---7-'' -"-J 

represents the discrete gradient associated to the chain with the Neumann boundary conditions. Its trans- 
pose operator is denoted (V'^)'^ : K^^^ — > M^. The particles are represented by their re-scaled positions 
q — (gi, g^v), so that the typical position and deviation of q is formally of order 1 with rrespectto N. 
This can be seen by considering particles in the chain as indexed by a; = e [0, 1]. The thermodynamic 
limit ^ oo of the Hamiltonian (j5.2p then converges to 



lim —HNiq,p)= / 7:p(a;) dx + ^^int 3- dx + Vcxt{qix)) dx , (5.4) 
N^oo N Jq 2 Jq \dxj Jq 

when X 1— > q{x) is a smooth function on [0, 1] satisfying the Neumann boundary conditions q'(0) — ^'(l) = 0. 
We choose to work with the scaling given by ()5.2p - ()5.4|) in order to prescribe the macroscopic timescale 
of the chain profile at order 1 with respect to N. Using the matrix notation, we write the stochastically 
perturbed equations of motion 

q = P 

fV7d\T, 



p = -(V'^)^<t(V'^q)-<,t(g)-7P + VA^^W^ 



with the fluctuation/dissipation identity acr^ — 2(3^^ j. 

Next we turn to the mass-penalized Hamiltonian model for this particle chain. The number of inde- 
pendent fast degrees of freedom that will be penalized is equal to — 1. The penalized degrees of freedom 
are given by the inter-particle distances S,i{q) = z.^ ^ qi — qi-i- The mass-penalization intensity has to scale 
formally as 

V := Vfq = vN 

to give a discrete operator meaning to the perturbed mass-matrix in 

Af,„ =Id -Ki>2(V'^)^M,V'^. 
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When AIz = Id the matrix perturbation is given by the discrete Laplacian 

Ad := -(V'')^V''. 

with the Neumann boundary conditions. FoUowing our general construction we obtain the mass-penahzed 
Hamiltonian 

^ N N 

i=l i=l 
and the stochasticaUy perturbed equations of motion with — Id are 

q = (Id - D^^a)-' (-(V'^)^<t(V''g) - <,t(<?) - 79 + CJ^/Nw) . (5.6) 

The form (|5.6p of the equations of motion clearly demonstrates the regularizing effect of the mass-matrix 
penalization. Indeed, all the high freciuencies generated by the microscopic forces —{y^Y'v[^^{y'^q) are 
filtered by a "low-pass filter" represented by the regularizing operator (id — i^^A^) ^ . 

5.2. Numerical results. In this section we demonstrate behavior of the IMMP method with the 
Hamiltonian (j5.5p for different values of the re-scaled mass-matrix penalization intensity v, the time step 
parameter (5t, and the system size N . In particular, we perform a systematic comparison with the exact 
dynamics (|5.2p integrated by the Leapfrog/ Verlet scheme. 

The simulations were carried out with the re-scaled quadratic exterior potential 

and with a double-well repulsive re-scaled interaction potential 

Vint{r) = (50 ((r - 0.1)^ - O.OS^))^ , r < 0.1 
Wint(r)=0, r>0.1, 

and the rescaled inverse temperature /3jv = 10 A^^"'^. We choose the penalizing mass-matrix to be the 
identity matrix Mz = Id , and the dissipation parameter 7 — 0.1. 

Test I: Comparison of macroscopic dynamics. 

We demonstrate approximation of the macroscopic dynamics by the IMMP method simulating two 
observables: (i) the length of the chain: I = qjf — qi, and (ii) the center of mass: c = qjv/2- The dynamical 
behavior for short times of the chain length I and the center of mass c for different values of the mass- 
matrix penalization intensity P is depicted in Figure 15.11 The same realization of the noise is used in each 
simulation, and no Metropolis rejection is used. One can observe the smoothing induced by the IMMP 
method, however, typical evolution timescales are not changed by the IMMP method. 

Next, the convergence to equilibrium of I and c are analyzed. The initial condition is taken out of 
equilibrium with zero potential energy and particles concentrated near the central position q = 0.5. The 
equilibrium probability density functions (PDF) of I and c are computed for the IMMP case and the exact 
dynamics in Figures l5.2i exact sampling ensuring identity of the two calculations. The PDFs are estimated 
using a Gaussian kernel estimator. The evolution of the relative entropy between a trajectory sample 
and the equilibrium distribution is plotted in Figure \5M We observe that the convergence of the IMMP 
dynamics is slightly faster than the exact dynamics with the Verlet integrator. Convergence to equilibrium 
on shorter times can be also observed by inspecting the auto-correlations in time series of l{t) and c(t) in 
Figure 15.31 Convergence of the IMMP method is again faster. However, it is interesting to note that the 
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mass ppenalization introduces an inertial artifact that produces oscillations in the auto-correlations of the 
center of mass evolution. 

The average transition time r(i/, N) of the chain center of mass c between q — 0.4 and q — 0.6 is 
summarized in Table [STTf for different values of the chain size N and the mass penalization P. The results 
show that the average time t{D, N) of the macroscopic dynamics of c scales as 0(1) with respect to the 
system size, and is continuous with respect to the exact dynamics when — > 0. 




Figure 5.1. Short time evolution of the chain length and the center of mass for different values of the IMMP mass- 
penalization p2 = {10-3, 10-*} and = {5.10-^,5.10-*}, N = 500. 





D'^ = IQ-'^ 


9'^ ^ 10-=* 


9^ = 10-* 


9'^ = IQ-'^ 


Leapfrog/ Verlet 


N = 


100 


1.6 (1) 


1.3 (1) 


1.2 (1) 


1.0 (1) 


1.0 (1) 


N = 


300 


1.5 (1) 


1.3 (1) 


1.2 (1) 


1.1 (1) 


1.0 (1) 


N = 


500 


1.4 (1) 


1.2 (1) 


1.2 (1) 


1.1 (1) 


1.0 (1) 



Table 5.1 



Average transition time of the chain center of mass from the position gjv/2 = 0-4 to the position gjv/2 = 0-6, simulated 
by the IMMP method. 



Test II: Relaxation of time-step restrictions. 



We first observe in Figure 5.5(a) the relationship between the time-step St and the mass penalty v. 



The average Metropolis acceptance rate is computed for different time-steps and different values of the 
mass-matrix penalization. The critical time-step, for which the average acceptance rate is approximately 
0.5, is significantly increased by the mass-matrix penalization, and is broadly proportional to the re-scaled 
mass-matrix penalty 9. 



In Figure 5.5(b) the critical time-step is studied for different values of the chain size N and of the mass- 
matrix penalization, and it is compared to the exact dynamics. Using linear regression the dependence 
of the critical time-step on the system size is estimated and compared to the exact Verlet integration in 
Table [5TW . These results are in a good agreement with Proposition 16 .31 that characterizes the dependence 
for the linear case. 



Conclusions. 



^Digits inside the brackets indicate the typical error of the last digit. 
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Figure 5.2. Equilibrium PDF of the chain length and the center of mass. Simulated by the IMMP method with the 
penalty = 10~^ . 



® IMMP =0.01 
♦ Leapt rog/VerIf 



© IMMP -0.01 
♦ Leapfrog/Verlel 



Figure 5.3. Convergence of the relative entropy of the PDF of a trajectory of the chain length and the center of mass 
towards equilibrium. Simulated by the IMMP method with the penalty = 10~^. The error bar represents extreme quartiles 
as computed over 200 realizations. 



The presented numerical studies demonstrate that for a prescribed macroscopic timescale of the chain, 
the IMMP aUows for arbitrary large (of order N) relaxation of the time-step restriction. This observation 
shows that for a large system, iV — > oo, the IMMP method can achieve arbitrarily large computational 
savings when compared to the standard Verlet algorithm used for integrating the exact dynamics. 

6. Numerical analysis of the particle chain. In this section, the particle chain model (|5.2p is 
analyzed in the special case of harmonic interactions. We consider the thermodynamic limit N +00 
where N is the size of the system. It is shown that the macroscopic dynamics of the IMMP method behaves 
continuously (uniformly with A'^, and in the L2 norm for the position profile) with respect to the re-scaled 
mass penalty parameter D. 

At the same time, the time-step stability of the IMMP numerical scheme (|4.1|) is compared with the 
standard Verlet scheme, and the critical time step is shown to be increased by a factor i?N . 

From the spectral point of view, the IMMP method behaves in this linear case as a low-pass filter. This 
proves, in this simplified case, the ability of IMMP method to respect macroscopic dynamical equivalence, 
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IMMP^O.OI 
Leapfrog/Verlei 




0.012- 
0.010- 



0.002- 
0.000- 
-0.00^ 
-0.00-1- 



IMMP -0.01 

Leapfrog/Verlel 




Figure 5.4. Auto- correlation of a trajectory of the chain length and the center of mass. Simulated by the IMMP method 
with the penalty 9^ = 10"'^. 



0.9- 
0.8- 
0.7- 
0.6- 
0.5- 



+ LeaplrogVerIf 

X IMMP^O.01 

® IMMP -0.1 

♦ IMMP-1 



+ IMMP -0.1 
X IMMP -0.01 
© LeapfrogVerlet 



(a) Acceptance rate 



(b) Critical time-step relaxation 



Figure 5.5. (a) The average acceptance rate of the Metropolis step for different values of the mass-matrix penalization 
= {1, 10"-*^, 10~^} and the system size N = 100. (b) Critical time-step that achieves the Metropolis acceptance rate of 0.5. 
Comparison between the IMMP dynamics (u^ = {lO""*^, 10~^}J, and the exact dynamics, for different values of the system 
size N. 



while saving computational time up to a factor of order 0{N). 

6.1. Conservation of macroscopic dynamics. We consider the penalized Hamiltonian (|5.5p when 
the interaction potential is quadratic Viat{f) = The penalizing matrix is the identity matrix Mz = Id , 

hence the penalized mass-tensor becomes = Id — P^Ad, and the fluctuation/dissipation tensor is taken 
proportional to the identity matrix. The system of equations ()5.6|) then becomes 



(6.1) 



(Id - p2Arf)-V 
Adq-v'^^,{q)-jq + aVNW, 



with fluctuation/dissipation identity 7 — 2a'^f3 ^. The associated canonical equilibrium distribution is 
then given by the re-scaled inverse temperature Pn = I3N~^. 
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z/"^ = 10 ^ = \{) Leapfrog/Verlet 


St cx iV-" 


a = 0.2(1) a = 0.2(1) a = 1.2(1) 



Table 5.2 

Scaling of the critical time-step for the Metropolis acceptance rate r = 0.5. 



In order to treat the limit N ^ oo, we introduce the £2-iiorm in the position space 

N 



ii2 1 2 It 



as well as the ft,_i-norm in the momentum space 



IbllL, 



N 2 , 

1=1 \ 4 = 1 



In the above expression, X]i=i Pi '^^^ seen as the orthogonal projection in on the one dimensional 
kernel of the Neumann discrete Laplacian A,;. H'zllo + llpll'ii is a quadratic form that endows the phase-space 
with a Hilbert space structure. 

Proposition 6.1 (Convergence of the macroscopic dynamics). Assume that the exterior potential 
Vcxt is bounded and that its derivative satisfies the Lipschitz condition 

Ikoxtfe) - ^^oxt(9i)IL_i < Ly \\q2 - qiWi^ , 

where Ly is independent of N . For any T > 0, let t {p'^ {t) , q'^ (t)) be the solution for t G [0,T] of the 
evolution equation i6.1]) with the initial condition 

(/(0),/(0)) = (AC;/V(0),g°(0)), 

where (p'^ (0) , q'^ (0)) is distributed according to the original equilibrium canonical distribution (associated 
with (|5.2p - (j5.3p ). Then for all t € [0, T] one has the uniform convergence 



lim lim sup I 



q-'{t)^q'^=^{t) 



0. 



Proof. We write X = ((7,p), and introduce the norm: 



\X\\.-. 



h-i 



The system (16. ip becomes a stochastic differential equation in the form 

dXf = ApXf + F(Xf ) + T.dWt , 

where by definition 



(6.2) 



A, := 



(Id - D^Ad)-^ 
Ad 



Fix) 







S := 





Na 



'v'^Al) ~ IP J 

Duhamel formula gives an implicit expression for differences of solutions of (|6.2p with the same noise 



Xf-X," = (e^^ 



{F{X^^)ds + Y.dWs) 



+e'^^{X^ - X"o) + f e^^(*-^)(i^(Xf ) - F{X^,))ds . 
Jo 



(6.3) 
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We estimate the individual terms on the right hand side in (|6.3p . We define P as the coordinate transfor- 
mation associated with the orthonormal spectral decomposition 

-Arf = p-Miag(5o,...,'5jv-i)F, 
where PP^ = Id , and the eigenvalues of the discrete Neumann Laplacian are given, for fc = 0, . . . , — 1, 

by 



AN^ sin^ 



kn \ 

2N J N^oo' 



2„2 



(6.4) 



Denoting the spectral coordinates 



we have 



X = {q,p) ^ {N-^''^Pp,N~^''^Pq) 



N-l 



4 -2 



N-1 



The spectral decomposition leads to a block diagonal form of the operator e^"* with diagonal 2x2 blocks 
in the spectral basis 



1 r 

1/ ' 



as well as for fc = 1 , . . . , — 1 







^(fc) ^(k) 
where is the 2x2 block associated with the coordinates {qk,Pk)- Since conserves the fc-mode 

energy dk{qk)'^ + (1 + i^^'^fc) "'^(Pfc)^, one can check that for any > 1 



||gA,t|||2 < 2 + 2^2. 



— 1/2 

Similarly, since in the sense of symmetric matrices Aft/^ < Id , we have the bound 



\F{Xn-F{X")l< M~y^vWin-vUq')+l{p'-p')) ^ <{Lp+^)\\X-^-X^^ 



Using random independence of Brownian increments we compute 



(^^ArAt-s) _ ^Ao(t-s)\^ ^^^^ 



t N 



ds . 



Applying Gronwall lemma in (|6.3p and collecting all terms we obtain 



E 



\x^-x^\ 



< 



X^-XX + mr 



where Ct is independent of N, and with being distributed canonically uit is given 



(6.5) 



N 



rriT = sup 
te[o,T] 
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For a given random vector X such that E 
imply 



\X\ 



< +00, Parseval identity and the inequahty 



Af-l 



E 



10 



^Apt „Aot\ 



X 



k=l 



2 



N-1 

<2^E 

fc=i 



X 



k,Q 



where ^ is the restriction to the k-th mode {qk,Pk)- Then one has, by orthogonahty of P, 



I < 



(6.6) 



N 







2 












*;,0 



^ 1 2 



Now the canonical distribution of X'^ — (9*^,^°) can be described as follows. Up to normalization, the 

distribution of X^ has the density e~''v ^»=i "'=''''^*» with respect to the Gaussian distribution with the 
covariance matrix /3~^Id for momenta variables, and the covariance matrix (/JA^)"^ for positions. Thus 
we have the bound 



as well as 



lim E 

N^oo 



E 



F{XO) 



XO 



k,0 



1 



k,0 



k,0 



where J-^ denotes the Fourier series expansion on [0, 1] with Neumann conditions, and {G^k)k>i are canonical 

1 

k^TT- 

the series is bounded 



centered Gaussian i.i.d. variables with the covariance matrix /? 



By the Lipschitz assumption 



k=l 



k.O 



< (Lv + 7)E 







2" 












o_ 



+00 

E 

fe=i 



2(L^+7) 



Since lim^^o |||e'^" * — *||| = 0, one can take the limit N — > +00 and use the uniform convergence of 
the series in (I6.6P to obtain limp^o l™Ar^+oo "t-t = in (j6.5p . The convergence of the initial condition 
limp^o liniAT^+oo E \\Xq — ^0 ||p follows by using similar arguments. The proof is complete. □ 

6.2. Relaxation of time-step stability restriction. To demonstrate improved stability properties 
of time integration algorithms we consider the IMMP scheme (14. 2|) associated with the mass-matrix penal- 
ized Hamiltonian ()5.5|1 . Note that when the constraints are linear, the leapfrog scheme (RATTLE) applied 
to an implicit Hamiltonian is identical to the usual leapfrog scheme for the associated explicit Hamiltonian 
(|5.5|) . We restrict the rigorous analysis to the quadratic interaction potential (fint(2/) = zero exterior 
potential (fext = 0), and to the mass-matrix penalization operator (Id — v^A^). The leapfrog scheme is 
defined as 

Pn+1/2 =Pn + — {-Ad)qn 
Qn+l = Qn + St 

Pn+1 =P„+l/2 + y(-Arf)g„+i . 
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Denoting the spectral variables for fc = 1, 



,N- 



1 + D^Sk 
Sk 



1/2 



1/2 



NPp 
NPq, 



(6.7) 



we write 



where 



Since det (Lfe) 



-'n+l 




-hk + Y 



and hk 



5t- 



(1 + ^^24) 



1/2 



1, the standard CFL stability condition is equivalent to 

|Tr(Lfe)| <2 

1. Thus we arrive at the following bound on the time 



which is fulfilled if and only if < 2 for all A: < 
step 



5t <2 mill 

0<A;<A' 



1 



1/2 



Summarizing the above calculations and recalling (|6.4p we have the following characterization of the 
stability properties. 

Proposition 6.2. Suppose Vcxt ~ and consider a harmonic interaction potential Wint(?') — ^ with 
the mass-matrix penalization M^^ = Id — i^^A^;. The leapfrog /Verlet integration of the IMMP h 



Hamiltonian (|5.2[) is stable in the spectral sense if and only if 

I 



1/2 



it < 



\ 



sin^ 



{N - l)7r 
27V 



armomc 



(6.8) 



Since we work with a Metropolis correction of the hybrid Monte-Carlo type, we are also interested in 
the limiting behavior of the energy variation compared to the temperature, i.e., 

PN{H{Pn+l,qn+l) - H{pn,qn)) , 

when {pmin) are distributed according to the canonical distribution. This quantity gives the average 
acceptance rate of the Metropolis correction. The result we present here is similar to [4j where the authors 
analyze infinite dimensional sampling with the standard Metropolis-Hastings Markov chains. 

Proposition 6.3. Suppose Vext = and consider a harmonic interaction potential Vint{r) = ^ with the 



mass-matrix penalization Mi, 



Id — Arf. Suppose the state variable X = q) is a random variable 



distributed according to the canonical distribution associated with the mass-matrix penalized Hamiltonian 
(15. 2p . Then the energy variation (Sn^H after one step of the leapfrog integration scheme converges in 
distribution, up to normalization and centering, to the Gaussian random variable 



jSN^H-mpf La 



■AA(0,1) 
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with the mean and variance in the infinite size asymptotics for the IMMP method v > and St — Stpf = o(l) 

N6t% ^ 2 NSt% 

and for the Verlet integration of exact dynamics with St = St^ = o{l/N) 



5 5 
niN ^ -N'^St%, and a% - -N"^ St% 

N~> + oc 8 N^+oc 4 



Proof. We start with a canonically distributed state X = {q,p), which is, by assumption on the form 
of the interaction potential, a Gaussian random vector. After changing to the spectral coordinates (|6.7p 
we have the spectral representation of the Hamiltonian 

N-l rl/2 

and introducing Gaussian random vectors U and V with the identity covariance matrix we can write 
X ^13 ' -jyj Uk, and w =/3 ' -j-rj Vk . 



We then compute explicitly the change of the Hamiltonian after one step of the leapfrog integration 



) ^']- (6-9) 



Since det (L^L^) = 1 the matrix Lj^L^ — Id has two positive eigenvalues (A^ — 1, 1/Afc — 1) which satisfy 

hi 



Afc + 1/Afe - 2 = Tr {LlLk ^ M ) = 

hl^ hi 



k - > -^g . 
(Afc - 1)2 + (1/Afe - 1)2 = Tr {LlLkf - 2Tr (L^L^) = + '-^ 



256 8 

Combing with (j6.9p we find 



12 



^-1 Ifi ^-1 /j6 ^ 

fe=i fe=i 

Moreover, the Lindenberg or simply Lyapunov condition in the general central limit theorem (see [H 
verified since we work with a sum of random variables, thus concluding the first part of the proof. 
Recalling 



h,.. = St-^-^^ 



k 7T \ 

N 2 ) 

'(4^ + ^^sin2(A|))i/2' 
we compute the convergent Riemann sums for p — 6 and p — 12. For the case i? we have 



stp 



N^oo iV ^ DP 
k=l 



lim — /i? — 
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If = we obtain 

N-l „,„ N-1 



lim -i— y hi = lim 5 y 2f sinP f 

fe=l fc=l ^ ^ 



„1 

= (5tP2P y sinP (|a;) c^a;- 

Thus for p = 6 we have that the series sums to 20St^. Then the asymptotic behavior foUows from the 
assumption St]^ <^ St%, and similarly 6t]^N^'^ < St%N^ in the case D — 0. U 

Remark 6.4. The two propositions proved in this section characterize the restrictions imposed by 
the stability of the resulting scheme. In Proposition 16. 2( stability in the large system size limit, N +oo, 
is equivalent to the inequality (|6.8p . In this case the restriction of the time-step size is imposed by the 
numerical integrator. On the other hand the stability for the scheme which uses a Metropolis corrector is 
linked to the acceptance rate of the Metropolis step. In Proposition 16.21 stability in the large system size 
limit is equivalent to the non-vanishing Metropolis acceptance rate, which is equivalent to bounded from 
above average energy variation mpf and bounded variance apf of the energy variation. In either case, the 
IMMP method {P > 0) induces a relative increase of order N for the boundary of numerical stability as 
compared to the exact dynamics 9 = integrated with the Verlet scheme. 

7. The stifT limit. Throughout this section, one introduces a potential function with an explicit 
dependence with respect to the fast variables {q, z) t—^ U (q, z) together with a stiffness parameter e. The 
potential energy V can then be written in the form 

Viq) = UiqM), 
e 

with a confining assumption 

inf U{q,z) > K{z), 

qeR'' 

where lim2_^oo K{z) = +oo strictly faster than at any logarithmic rate. The functions ^ are then indeed 
"fast" degrees of freedom (fDOFs) in the limit e ^ 0, the system being confined to the slow sub-manifold 

Mo^{q\a<l)^0}- 

We will prove, that under appropriate scaling of the mass penalty i^^ = the IMMP method is 
asymptotically stable in the stiff limit, converging towards standard effective dynamics on the slow manifold 
Mo. 

7.1. Thermostatted stiff systems. The canonical distribution becomes 

f,,{dpdq) = ^e-'^(^p"M-V+t^(?,^))dpdg. (7.1) 

In the infinite stiffness limit (e 0) the measure concentrates on the slow manifold A4o- The limit is 
computed using the co-area formula (see Appendix [A] for relevant definitions of surface measures). In order 
to characterize the limiting measure we introduce the effective potential 

UM - e-'3^(''") dz . (7.2) 

Lemma 7.1. In the infinite stiffness limit (e — > OJ, the highly oscillatory canonical distribution (j7.ip 
converges /ie^/^o ftn distribution) towards ^.^{dpdq), which is supported on M.q, and defined as 

Mdpdq) = ^e-/3(^p"M-V+^4«(<?)) dpS^^,)=o{dq) . (7.3) 
^0 
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Its marginal distribution in position is given, up to the normalization, by 

e-Py^«('^)6^^,^Mdq). (7.4) 

Proof. It is sufBcient to consider distributions in the position variable q only. Let be a S- 
neighborhood of A^o where 

dq = e"(55(,)=„(d(7) dz . 

We construct a decomposition (p = (pi + ip2 oi continuous bounded observables such that supp pi C 14^ 
and supp(p2 H U^^'^ — 0. Using the confining property of U{q, •) we obtain 

^(q)e'PUi,,^)^q = 1 (/.i(<z)e-'3^(«-)<5^(,).,,(dg) dz + 0{e-^^<^'^'^^) . 
By continuity of e i-^- J </'i('?)6~''^^*'^^'5^(g)=ez('^'?) ^^^id by the dominated convergence theorem 

VPi(9)e-'^^('''^)<5^(,)=,,(dg)dz^ J ^i(g)e-^^-(«)<5^(,)=o(d<z) = J v{q)e-f'''--^^H^(^)=,{dq) , 

and the result follows after normalization. □ 

The infinite stiffness limit (e 0) of highly oscillatory dynamics has been studied in a series of papers 
[551 [571 [771 [51 [5^ [53] . The limiting dynamics can be fully characterized in special cases. For example, when 
the highly oscillatory potential is linear and non-resonant (at least almost everywhere on the trajectory, 
see [37]), it can be described through adiabatic effective potentials. See also [lOl [6] for a recent work on 
some related numerical issues. However, when the system is thermostatted, one can postulate an "ad hoc" 
effective dynamics ([55]) exhibiting the appropriate limiting canonical distribution given by (|7.4p . Such 
dynamics can be obtained by constraining the system to the slow manifold M.o, and adding a correcting 
entropic potential, sometimes called "Fixman" corrector from [191, which is due to the geometry of A^O; 
and is given by 

^fix(g) = :^ln(detG(<z)) , (7.5) 



where G{q) is the n x n Gram matrix defined in 

In general, since the effective potential (17. 2p is not explicit, one may need to couple the system with 
virtual fast degrees of freedom to enforce the appropriate effective dynamics associated with (17.21) . The 
resulting extended Hamiltonian is then defined on the state space T* {Aio x R") (the cotangent bundle) 
and is given by 

H,ff{p,p„ q, z) = ip^Af- V + ipjM^V. + U{q, z) + yfi,(g) 

e(<7)-0. {C) 

Definition 7.2 (Effective Langevin process with constraints). The constrained Langevin process 
associated with Hamiltonian (|7.6p is defined by the following stochastic differential equations 

' q = M^V 
z = M-V. 

p - -ViU{q, z) - VqVfiM -iq + dW- V2i A (7.7) 
Pz = -VzU{q,z) - + a^W^ 
1^(9) = 0, (C) 
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where Vi and V2 are respectively derivatives with respect to the first and second variable of the function 
U{q,z), W (resp. ) is the standard multi- dimensional white noise, 7 (resp. jz) adx d (resp. n x n) 
symmetric positive semi-definite dissipation matrix, a (resp. a^) is the fluctuation matrix satisfying aa^ = 
^7 (resp. (TzCrJ = '^Iz)- The processes A € R" are Lagrange multipliers associated with the constraints 
(C) and adapted with respect to the white noise. 

We formulate reversibility of this process as a separate lemma. 

Lemma 7.3. The process defined in |7. ?[ ) is reversible with respect to the associated canonical distri- 
bution whose marginal distribution in {q,p) variables is 

l,,s{dpdq) = ^e-''(^f"^^"'^+^-(«)+^-(«))aT.A.o(dpdg) (7.8) 

^eff 

with the q-marginal 

e-Py^«i<^) 5^(,)=,{dq) . 

When 7 and 72 are strictly positive definite, the process is ergodic. 

Proof. The process (|7.7p is a Langevin process with mechanical constraints, exhibiting reversibility 
properties with respect to the associated Boltzmann canonical measure (see the summary in Appendix |BJ . 
Then the q-marginal is obtained by remarking that the integration of any function of ^p^ M~^p-\-^p'^M~^pz 
with respect to dpz cyT^Mo{dp) results in a constant independent of q. □ 

The properties of thermostatted highly oscillatory systems are summarized in Table W7\\ 





Finite stiffness 


Infinite stiffness limit 


Infinite stiffness 




e > 





e = 


Dynamics 


Highly 


Adiabatic 


Effective with 




oscillatory 


(if non-resonant) 


constraints 




+ fluct./diss. 


+ non-Markov fluct./diss. 


+ ffuct./diss. 


Statistics 


Canonical 


Positions on A^o, 
free velocities. 


Canonical on T*Mq, 
geometric corrector. 


Numerics 


Leap frog/ Verlet 


Time-step 


Leapfrog/Verlet with 




+ fluct./diss. 


restrictions {5t = o(e)) 


constraints 
-1- ffuct./diss. 



Table 7.1 



Stiff Hamiltonian systems and associated commonly used numerical methods (A4o denotes the slow manifold). Two 
different schemes are required for the stiff system and its effective Markovian approximation. 

7.2. Stability of the IMMP dynamics. We assume that the mass-matrix penalty parameter v = 
grows to infinity in such a way that 

lim ej/f = D . 

The original Hamiltonian with the stiffness parameter is expressed explicitly as 

He{p. q) = \p^M-^p + U{q, ^) , (7.9) 
and including the mass-matrix penalization one gets 

H.AP.,,q) = \pIM-^p,^ -f U{q, ^) + l^fi.,.,(g) , (7.10) 
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or in its implicit formulation 

11 z 

2 2 ..e ^^^^^^ 

m = -z- (C.J 

One immediately sees that iJiMMP is non-singular when e — s- and converges to the effective Hamiltonian 
on the slow manifold, 

i/eff,p(g, z,p,p,) = ^p'^M-^p + ^p'^M-^p, + U{q, |) + Vfi^iq) 

The expression (j7.1ip represents a minor generalization of H^s (j7.6p . but it leads to the same canonical 
marginal distribution ^osidpdq) in (j>, q) variables as given by (|7.8p . The continuity in e of -ffiMMP implies 
stability of the associated dynamics and their numerical integrators. We first derive the limits of the 
original and penalized canonical distribution. 

Proposition 7.4 (Limits of canonical distributions). Consider the canonical distributions fj,^^{dpdq 
associated with the mass penalized Hamiltonian (j7.10p . but considered with respect to the (jp = MM~^p^^, q) 
variables. In the sense of weak convergence of measures we have fi^^ ^/^off as e — > with ficft defined by 

(EHl). 

Proof. The first convergence is proved in Lemma 1 7. II For the second one, the following notation will 
be used 

J Sg,ez{dq) ^ S^(^q)^^^{dq) 

To prove the convergence towards /ioff we consider a (5-neighborhood lA^ of A^o where 

dpdq = jg^j^^ ^q^ez{dq) dz Sp^ep, (dp) dpz , 

and a decomposition of the bounded observable (in {p, q) variables) ip — ipi + ip2 such that supp ipi C lA^ 
and supp</?2 nZ^"^/^ = 0. Thus, keeping in mind that p^^ — M^^AI^^p, and using the confining property of 
the potential U{q, ■) we obtain 

J ^{p,q)e-P"-'^dp,, dq^ J vJi(p,g)e-^^-rfp,, dg + 0(e-'^^(*/'=)) = + ©(e-''^^^/^^)) . (7.13) 

Applying the change of variables p,^^ — Mi^^M^^p yields 

dp^^ = det (M^^ M-^) dp = vf' det det (G + ^M"^) dp , 

and setting eM~^pz = V M^^p and ez = ^{q) we get 

1 z 

HyXVi^e-,q) = 7:P^ M^^p + v'^e'^p^ M^^p^_ + U{q, — ) + Vax.i/.(g) = HiuMpiq, z,p,Pz) • 

2 l^ei 

Thus substituting back to (|7.13p we obtain 

= (i/,e)2« / ^ie-'5^™"-(«-^-P-f^)det(G+^M-i)(5p,,p^(dp)dp,<5,,e.(dg)dz, 
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and thus 




ifie 



/3ff.„,.(9,.,p,P.)^g^ (G) (dp) dp, 5a^)^^{dq) dz . 



Using the co-area formula we obtain 



det (G) 5xj(^(q)M-^p=a{dp)5(^(^q)^a{dq) = <7T'Ma{dpdq) , 



which leads to the final result after integration of the (p^, z) variables and normalization. □ 

Remark 7.5. Due to the fast oscillations, the distribution of impulses in the limiting distribution 
/iQ in (|7.3p is uncorrelated, whereas after the mass-matrix penalization, the limiting distribution (j3.5p has 
almost surely co-tangent impulses (i.e., satisfying the constraints Vq^M^^p — 0). This explains the role 
of the corrected potential energy Vrx taking into account the curvature of A^o- 

In the next step we inspect the infinite stiffness asymptotic of the penalized dynamics. 
Proposition 7.6 (Infinite stiffness limit). When e with v = ^ ^ and V{q,£_{q)) = U{q, ^), 
the IMMP Langevin stochastic process p.6p converges weakly towards the following coupled limiting pro- 
cesses with constraints 



where Vi and V2 are respectively derivatives with respect to the first and second variable of the function 
U{q,z), and {At}t>o are adapted stochastic processes defining the Lagrange multipliers associated with the 
constraints (C). 

The process {(?t,Pt}t>o defines an effective dynamics with constraints (Definition \7. 2^ for thermostatted 
highly oscillatory systems. It is reversible with respect to its stationary canonical distribution given by figS 
(j7.8p . and is ergodic when (7,72) are strictly positive definite. 

Proof. The proof is similar to the proof of Proposition 13.71 Here we have 



< 



'q = M-'p, 

p = -WiU{q, 4) - %Ffi,(g) ^^q + aW- %CA , 
1/ 

^(9) = 0, 

1 z 

Pz = --V2U{q, -) - JzZ + <JzWz . 



(G) 



(7.14) 



VqU = VlU + -Vq^V^U , 



and (j3.6[) translates, up to a change of Lagrange multipliers, into 



' q = M-^p 
i = M^^Pz 



P 



Vit/ - \/qVti^.^{q) --fq + aW- V,e A 



(7.15) 



< 



Pz 



— -jzZ + azWz + — 



Z 



(G.J. 



The rest follows the proof of Proposition 13. 71 
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□ 

Remark 7.7. When v +00, by a classical averaging argument (see, e.g., [27]), one can check that 
the limiting dynamics are the effective dynamics pointed out in [33] 



' q = Ar^ 



P 



.e(9) = o. 



(7.16) 



(C) 



with the stationary canonical distribution (|7.8 



7.3. Stability of the IMMP integrator. The numerical scheme (Scheme 14.11 proposed for the 
IMMP method p.6p ) is also stable in the limit of infinite stiffness e 0. Recall that we consider a 
reversible, measure preserving numerical flow ^gl{p,Pz, q, z) associated with Hamiltonian (|7.1ip iJiMMP 
with constraints (modified potentials could similarly be considered). 

Proposition 7.8 (Asymptotic stability). In the limit ev^ —f v, the numerical flow <i>^^ associated 
with the leapfrog /Verlet integrator with constraints for the IMMP Hamiltonian (j7.1ip converges towards 
the numerical flow <i>^j, which is the leapfrog /Verlet integrator with geometric constraints associated with 
effective Hamiltonian (|7.12p on the slow manifold. 

Proof. The statement is a direct consequence of the implicit function theorem and the continuity of 
the leapfrog integrator with constraints (|4.2|) with respect to the parameter D = ev^. Indeed, considering 
the shift of Lagrange multipliers 



A ^ A + -V2U 
e 



and taking the limit e ^ we obtain the appropriate leapfrog scheme 



Pn+1/2 =Pn- ^ViC/(g„, ^) - V,^(g„)A„+i/2 



2 

St. 



V 



Pn+1/2 =Pn- ^V2[/(g„, y) 

= 9n + StM^'^Pn+i/2 
Zn+l = Zn + 5tM-^pl^^/^ 

e('7«+i) = o 

Pn+l ^Pn+1/2 - ^Vif/((7„+i, -^) - Vq^(g„+i)A„+i 



Pn+1 — Pn+l 



|^V2C/(g„+i,^) 



V 



I V,e(<Z«+i)M-Vn+i =0. 



(Cl/2 



By convergence of the Hamiltonian (|7.1ip to (|7.12p , similar asymptotic stability properties holds when 
a Metropolis step is introduced. 

The results and properties discussed in this section are summarized in Table 17.21 
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Zero mass 


Positive 


Infinite 




penalization 


mass-penalization 


stiffness limit 






e,iy>0 




Dynamics 


Highly oscillatory 


IMMP 


Effective with 




+ fluct./diss. 


+ fluct./diss. 


constraints+ fluct./diss. 


Statistics 


Canonical 


Canonical with 
correlated velocities 


Canonical on 
T*Ma 


Numerics 


IMMP + fluct./diss. 



Table 7.2 



The IMMP dynamics and the Verlet numerical integration are both asymptotically stable in the infinite stiffness regime 
if J ^ 9 < +00. // the mass-penalization vanishes (v = 0) one recovers the original physical stiff system. The canonical 
distribution is always exact in the position variable. Notice that due to the penalized mass-matrix (v > 0) the statistics have 
correlated velocities. 



Appendix A. Surface measures. Let M be endowed with the scalar product given by the positive 
definite matrix Af , and consider M.^ a family of sub-manifolds of co-dimension n implicitly defined by n 
independent functions AAz — {x € M. | ^i('i') = zi, ..,^„((7) = z„} for z in a neighborhood of the origin. For 
each z in a neighborhood of the origin the conditional measure (5^(g)^^(c?(7) is a measure on AA^ defined in 
such a way that it satisfies the chain rule for conditional expectations with respect to the Lebesgue measure 
dq, i.e., 

dq ^ 5(^(q)=,Xdq) dz . (A.l) 

The surface measure o't^mA^p) Hausdorff measure induced by the metric on the co-tangent 

space T*M:, = {p I VjC(g)M-V = 0}; and in the same way, (JMAdl) is the Hausdorff measure induced 
by the metric M on the sub-manifold A4z- It is important to note that, although this is not explicit in 
notation, a is defined with respect to the mass-tensor M of the mechanical system. The Liouville measure 
CTT'M^idpdq) on the co-tangent bundle T*A4z is the volume form induced on 

T*M. = {{p, q) I v^e('z)M-V = 0, aq) = A 

by the usual symplectic form dp A dq. It can be described in terms of surface measures as follows 

c^T-M, id-pdq) = (TTjM, id-p) {dq) ■ 

Finally, the co-area formula (see [17] for a general reference) defines the relative probability density 
between (5^(q)^^(o?q) and aMAd<l)- 

Proposition A.l (Co-area formula). Given the invertible Gram matrix associated with the constraints 
^{q) — z in a neighborhood of = {ll^il) = z} 

G(<7) = V^eM^'V,C, 

one has 

S,,,,..idq) = -^^==aMAd<l) . 



Appendix B. Langevin processes. Defining the Poisson bracket 
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and the dissipation tensor 

d(g) = <Jq , 

where cr is the fluctuation matrix in Definition 12.11 the Markov generator of the Langcvin process in 
Definition 12.11 is 

C^{-,H} + ^{d,{d^,-}e-^"}e^". 

The generator £ satisfies 

j >PiC{^2)e-^" dpdq = J C*{ipi)ip2e-I^"dpdq, 

where 

The generator C* defines a Langevin process with the time-reversed Hamiltonian (—H). Reversibihty of the 
process implies that the canonical measure is stationary. Furthermore, if the initial state of the system is a 
canonically distributed random variable, the probability distribution of a trajectory after the time-reversal 
is given by a Langevin process with the generator C* . When H has the form H{p, q) = ^p^M~^p + V{q), 
reversal of impulses {p —^ —p) leads to time-reversed dynamics, and a process with generator C* can be 
constructed by the following simple steps: 

1. Reverse momenta (p — p). 

2. Draw a random path with generator C. 

3. Reverse again momenta {p —p)- 

When holonomic constraints, for instance, of the form 

are introduced, it is useful to define the Poisson bracket on the co-tangent bundle T*A4z 

a, 6 

where F is the symplectic Gram matrix of the full constraints 

As a basic result of symplectic geometry (see [T] ) , one recovers the divergence formula with respect to the 
bracket { ■ , ■ }^ and the Liouville measure (Tt'-m^ (dpdq) 

J {■ r}M,'^T-MAdpdq) = 0. 

Given a constrained Langevin process in a stochastic differential equation form 

q = VpH , 

p = -\/qH - 79 + rrVF - VgC A , 
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where A are Lagrange multipliers associated with the constraints ^(g) = 0, adapted with respect to the 
noise W, the process {pt,qt}t>o obeys hidden velocity constraints and is characterized by the stochastic 
differential equations 

p = -VqH --fq + aW - VqS A , 

where A are Lagrange multipliers associated with the full constraints q) = 0. The Markov generator 
of this process can be written in the form 



Cm.={-, H}^^ + i {d, {d^, • e-^"}^ e 



demonstrating the reversibility with respect to the constrained canonical measure e ^^aT*M^{dpdq). 

Appendix C. Exact sampling of fluctuation/dissipation perturbations. In this section, we 
recall how to perform exact sampling of fluctuation/dissipation perturbations. Since we only work with 
impulses, we refer to the system by using the impulse variables p only. Note that throughout the paper, we 
also use extended variables {p,pz), however, the presentation that follows covers general cases. The kinetic 
energy of the system is ^p^AIp. We impose constraints p-^M^^V^^ = on impulses, thus p G T*M and 
hence the associated orthogonal projector on T*A4 is 

P = Id - Wq^G-'^W^^M-^ . 

The stochastic differential equations of motion on impulses that are integrated on a time-step interval are 

(p = -jA4-'p + aW-Wg^X, 

with the usual fluctuation/dissipation relation aa'^ = 2/3~^7. The Gaussian distribution of impulses 

le-^p^''-'PaT^.Midp) (C.2) 

is invariant under the dynamics (|C.1|) . 

Proposition C.l (Exact sampling of stochastic perturbation). Given the mass matrix M, suppose 
either 6t or 7 are small enough so that the condition 

yM-i<7 (C.3) 

holds in the sense of symmetric semi-definite matrices. Let U be a centered and normalized Gaussian 
vector. Consider the mid-point Euler scheme with constraints 

Pn+l ^Pn - y7^~^(Pn + "/StaU - Vg^ A„+i 

where A„+i is the Lagrange multiplier associated with the constraint {Cp). The Markov kernel defined by 
the transition Pn Pn+i is reversible with respect to the Gaussian distribution (jC.2[) . 

Proof. After calculating the Lagrange multiplier the expression (|C.4p can be written as 

Pn+i = Pn - jP-rP^M-\p„+p,,+i) + VStPaU. 
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Consider the new variable p = ^^^p, and define the symmetric matrix 

L = yM~l/2p^pT^^-l/2 ^ 

as well as K, such that KK^ = L. In terms of these new variables we obtain from (jC.4p 

p„+i = (Id +L)-\ld -L)pn + 2{ld +L)-^KU. (C.5) 

Moreover, the product measure (^t* Mo{dpn) (Jt'A^o ('^Pn+i) t^*^ measure induced on the linear subspace 
of constraints by the scalar product M^^ and the Lebesgue measure dpn dpn+i- Thus in the variables 
(PmPn+i) this measure becomes, up to a constant, the measure induced by the usual Euclidean structure. 
As a consequence the log density of the random variable (pn,Pn+i) defined by (jC.5[) with respect to this 
latter measure is equal to 

-l\Pn\''-lipn+l-ild +Ly\ld -L)pnf L-\ld + i)^ (p„+i - (Id + L)-\ld -L)pn) 

Z o 

= -lpl+lL-\ld + Lfpn+l ~ \fnL-\ld + Lfpn , 

which is indeed symmetric between p„ and Pn+i- Hence we have shown the reversibility of the induced 
Markov kernel and consequently stationarity of the canonical Gaussian distribution. □ 
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